% replication_code_table_1_3.m
% 
% This matlab code replicates tables 1 and 3 in Miyamoto, W., T.L. Nguyen, 
% and H. Oh (2023), "In Search of Dominant Drivers of the Real Exchange 
% Rate," Review of Economics and Statistics.
%



%% setup
clear all;close all;clc
setup_workspace;


%% load files and add UIP wedge computation
jj_case=[814:820]; % baseline max share VAR cases stored for each country
% Case:    { 814,  815,  816,  817,  818,  819,  820 }
% Country: {'USA','CAN','JPN','GBR','DEU','FRA','ITA'}

jjj=0;
for jj=jj_case
    jjj=jjj+1;
    
    % set folder location to load VAR results
    ps_tmp.fpath.cucd=pwd; % current folder
    ps_tmp.file_name=['var_case' num2str(jj)];
    cd([ps.fpath.Case '/' ps_tmp.file_name]);
    load([ps_tmp.file_name '_var_result'])
    cd(ps_tmp.fpath.cucd)
    
    setup_workspace;
    
    % compute FEVD
    opt.T_fevd=30+1;%21; % fevd horizon length
    opt.var_fevd=[1:ps.N+1]; % variable index to compute fevd
    
    % Variable ordering in VAR:
    % 1 {'Output'} 2 {'Consumption'} 3 {'Hours Worked'} 4 {'NXY'} 
    % 5 {'RER'} 6 {'Inflation'} 7 {'Interest Rate'} 8 {'Investment'}
    % 
    % we add: 9 {'UIP wedge'}, using the current and future RER, inflation,
    % and interest rate dynamics
    
    % UIP wedge shock options:
    opt.F1=[ 0 0 0 0  1 -1 0 0];
    opt.F2=[ 0 0 0 0  -1 0 1 0]; 

    ii=0;
    % case of MBC identification : k-th variable
    for k_th=[1 2 3 4 5 6 7 8] %:ps.N % k-th variable (dominant shock variable)
        ii=ii+1;
        Y_fevd_draws=zeros(ps.N+1,opt.T_fevd,ps.nsave); % preallocation

        % variable by horizon by draw
        for i=1:ps.nsave % draws

            % companion form coefficient hx : X(t)=hx*X(t-1)+netashock*e(t)
            hx=zeros(ps.N*ps.p,ps.N*ps.p); % preallocation
            hx(ps.N+1:end,1:ps.N*(ps.p-1))=eye((ps.p-1)*ps.N);

            tmp=out.ALPHA_draws(2:end,:,i); % take out constants
            hx(1:ps.N,1:ps.N*ps.p)=transpose(tmp); % transpose 

            % create shock coefficient netashock 
            netashock=zeros(ps.N*ps.p,ps.N);
            tmp_chol=chol(out.SIGMA_draws(:,:,i),'lower');
            netashock(1:ps.N,1)=tmp_chol*rs.q_draws{k_th}(:,i); % first shock

            netashock_all=zeros(ps.N*ps.p,ps.N);
            netashock_all(1:ps.N,:)=tmp_chol;

            % set control variables 
            gx=eye(ps.N,ps.N*ps.p);
            % extend gx to include the UIP wedge variable
            In=zeros(ps.N,ps.N*ps.p);
            In(:,1:ps.N)=eye(ps.N);

            gx_ext=[gx; opt.F1*In*hx+opt.F2*In];

            % compute fevd
            Y=fevd_q(gx_ext,hx,netashock,netashock_all,opt);
            Y_fevd_draws(:,:,i)=Y.fevd;
        end
        
        % compute median fevd for each variable
        fevd_y(:,ii)     = median(Y_fevd_draws(1,:,:),3);
        fevd_c(:,ii)     = median(Y_fevd_draws(2,:,:),3);
        fevd_h(:,ii)     = median(Y_fevd_draws(3,:,:),3);
        fevd_nxy(:,ii)   = median(Y_fevd_draws(4,:,:),3);
        fevd_rer(:,ii)   = median(Y_fevd_draws(5,:,:),3);
        fevd_inf(:,ii)   = median(Y_fevd_draws(6,:,:),3);
        fevd_int(:,ii)   = median(Y_fevd_draws(7,:,:),3);
        fevd_inv(:,ii)   = median(Y_fevd_draws(8,:,:),3);
        fevd_wedge(:,ii) = median(Y_fevd_draws(9,:,:),3);
    end
    % for each country, store median fevd for each variable
    fevd_y_all(:,:,jjj)     = fevd_y;
    fevd_c_all(:,:,jjj)     = fevd_c;
    fevd_h_all(:,:,jjj)     = fevd_h;
    fevd_nxy_all(:,:,jjj)   = fevd_nxy;
    fevd_rer_all(:,:,jjj)   = fevd_rer;
    fevd_inf_all(:,:,jjj)   = fevd_inf;
    fevd_int_all(:,:,jjj)   = fevd_int;
    fevd_inv_all(:,:,jjj)   = fevd_inv;
    fevd_wedge_all(:,:,jjj) = fevd_wedge;
    
end


%% compute results for tables 1 and 3

% FEVD: h=4
tmp_h  = 100*squeeze(fevd_y_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_y_h4 = [tmp_h_med; tmp_h]; % output

tmp_h  = 100*squeeze(fevd_c_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_c_h4 = [tmp_h_med; tmp_h]; % consumption

tmp_h  = 100*squeeze(fevd_h_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_h_h4 = [tmp_h_med; tmp_h]; % hours worked

tmp_h  = 100*squeeze(fevd_nxy_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_nxy_h4 = [tmp_h_med; tmp_h]; % net exports to output ratio

tmp_h  = 100*squeeze(fevd_rer_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_rer_h4 = [tmp_h_med; tmp_h]; % real exchange rate

tmp_h  = 100*squeeze(fevd_inf_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_inf_h4 = [tmp_h_med; tmp_h]; % inflation

tmp_h  = 100*squeeze(fevd_int_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_int_h4 = [tmp_h_med; tmp_h]; % interest rate

tmp_h  = 100*squeeze(fevd_inv_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_inv_h4 = [tmp_h_med; tmp_h]; % investment

tmp_h  = 100*squeeze(fevd_wedge_all(1+4,:,:))'; tmp_h_med = median(tmp_h); 
fevd_wedge_h4 = [tmp_h_med; tmp_h]; % UIP wedge


% FEVD: h=20
tmp_h  = 100*squeeze(fevd_y_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_y_h20 = [tmp_h_med; tmp_h]; % output

tmp_h  = 100*squeeze(fevd_c_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_c_h20 = [tmp_h_med; tmp_h]; % consumption

tmp_h  = 100*squeeze(fevd_h_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_h_h20 = [tmp_h_med; tmp_h]; % hours worked

tmp_h  = 100*squeeze(fevd_nxy_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_nxy_h20 = [tmp_h_med; tmp_h]; % net exports to output ratio

tmp_h  = 100*squeeze(fevd_rer_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_rer_h20 = [tmp_h_med; tmp_h]; % real exchange rate

tmp_h  = 100*squeeze(fevd_inf_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_inf_h20 = [tmp_h_med; tmp_h]; % inflation

tmp_h  = 100*squeeze(fevd_int_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_int_h20 = [tmp_h_med; tmp_h]; % interest rate

tmp_h  = 100*squeeze(fevd_inv_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_inv_h20 = [tmp_h_med; tmp_h]; % investment

tmp_h  = 100*squeeze(fevd_wedge_all(1+20,:,:))'; tmp_h_med = median(tmp_h); 
fevd_wedge_h20 = [tmp_h_med; tmp_h]; % UIP wedge



% Create FEVD table
fevdtable_shock1_h4 = array2table([fevd_y_h4(:,1) fevd_c_h4(:,1) fevd_inv_h4(:,1) ...
    fevd_h_h4(:,1) fevd_nxy_h4(:,1) fevd_rer_h4(:,1) fevd_inf_h4(:,1) fevd_int_h4(:,1) fevd_wedge_h4(:,1)],...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate','UIP Wedge'},...
    'RowName',{'Median','USA','CAN','JPN','GBR','DEU','FRA','ITA'});

fevdtable_shock1_h20 = array2table([fevd_y_h20(:,1) fevd_c_h20(:,1) fevd_inv_h20(:,1) ...
    fevd_h_h20(:,1) fevd_nxy_h20(:,1) fevd_rer_h20(:,1) fevd_inf_h20(:,1) fevd_int_h20(:,1) fevd_wedge_h20(:,1)],...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate','UIP Wedge'},...
    'RowName',{'Median','USA','CAN','JPN','GBR','DEU','FRA','ITA'});

fevdtable_shock5_h4 = array2table([fevd_y_h4(:,5) fevd_c_h4(:,5) fevd_inv_h4(:,5) ...
    fevd_h_h4(:,5) fevd_nxy_h4(:,5) fevd_rer_h4(:,5) fevd_inf_h4(:,5) fevd_int_h4(:,5) fevd_wedge_h4(:,5)],...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate','UIP Wedge'},...
    'RowName',{'Median','USA','CAN','JPN','GBR','DEU','FRA','ITA'});

fevdtable_shock5_h20 = array2table([fevd_y_h20(:,5) fevd_c_h20(:,5) fevd_inv_h20(:,5) ...
    fevd_h_h20(:,5) fevd_nxy_h20(:,5) fevd_rer_h20(:,5) fevd_inf_h20(:,5) fevd_int_h20(:,5) fevd_wedge_h20(:,5)],...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate','UIP Wedge'},...
    'RowName',{'Median','USA','CAN','JPN','GBR','DEU','FRA','ITA'});

disp('-------')
disp('Table 1')
disp('-------')
disp('Dominant relative output shock; h=4')
disp(fevdtable_shock1_h4)
disp(' ')
disp('Dominant relative output shock; h=20')
disp(fevdtable_shock1_h20)
disp(' ')
disp('Dominant real exchange rate shock; h=4')
disp(fevdtable_shock5_h4)
disp(' ')
disp('Dominant real exchange rate shock; h=20')
disp(fevdtable_shock5_h20)
disp(' ')


% Create FEVD table
fevdtable_rer_h4 = array2table([fevd_rer_h4(:,1) fevd_rer_h4(:,2) fevd_rer_h4(:,8) ...
    fevd_rer_h4(:,3) fevd_rer_h4(:,4) fevd_rer_h4(:,6) fevd_rer_h4(:,7)],...
    'VariableNames',{'y shock','c shock','inv shock','h shock','nxy shock','inf shock','int shock'},...
    'RowName',{'Median','USA','CAN','JPN','GBR','DEU','FRA','ITA'});

fevdtable_rer_h20 = array2table([fevd_rer_h20(:,1) fevd_rer_h20(:,2) fevd_rer_h20(:,8) ...
    fevd_rer_h20(:,3) fevd_rer_h20(:,4) fevd_rer_h20(:,6) fevd_rer_h20(:,7)],...
    'VariableNames',{'y shock','c shock','inv shock','h shock','nxy shock','inf shock','int shock'},...
    'RowName',{'Median','USA','CAN','JPN','GBR','DEU','FRA','ITA'});


disp('-------')
disp('Table 3')
disp('-------')
disp('RER FEVD to each dominant shock; h=4')
disp(fevdtable_rer_h4)
disp(' ')
disp('RER FEVD to each dominant shock; h=20')
disp(fevdtable_rer_h20)
disp(' ')

